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Formulation  of  a  Fast  Iterative  Method 
for  Transonic  Flow  Calculations 


Reliable  but  slow  methods  for  calculating  transonic 
flows  have  been  developed  in  recent  years  [1,2,3].   These 
use  central  difference  formulas  in  the  subsonic  zone  and 
upwind  difference  formulas  in  the  supersonic  zone  to  ensure 
the  proper  region  of  dependence  and  jump  conditions.   The 
resulting  difference  equations  are  then  solved  by  an  itera- 
tion procedure  derived  from  the  method  of  successive 
overrelaxation.    This  note  describes  results  obtained  by 
using  a  fast  Poisson  solver  to   accelerate  the  rate  of 
convergence  of  the  iterative  scheme. 


Note:    This  work  was  also  partially  supported  by  NASA 
Grants  NGR-33016-167  and  -201. 
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It  was  proposed  by  Martin  and  Lomax  [4]  that  a  fast 
elliptic  solver  could  be  used  to  generate  an  iterative 
scheme  for  solving  the  difference  equations  appearing  in 
compressible  flow  calculations.   In  the  simplest  case 
consider  the  small  disturbance  equation 

(1-M^)(})    +4)    =0 
XX    yy 

where  (^      is  the  velocity  potential,  and   M   is  the  local 
Mach  number,   which  is  related  to  the  free  stream  Mach 
number   M   by  the  formula 

oo     -' 

M   =  M^  (1  +  (y+1)  (J)  )  . 

An  iterative  scheme  can  be  constructed  by  putting  the 
Laplacian  on  the  left  and  the  nonlinear  terms  on  the  right. 

Let   V   be  the  solution  for   *   at  the  nth   iteration.  Then 

n  ^ 

2 

Av  ,  ,  =  M   — t"  V 
n+1       „  2   n 
8x 

To  see  that  this  can  be  expected  to  converge  consider  the 

2  2 

linearized  equation  where   M    is  replaced  by   M  .   If  P 

and  Q   are  nonnegative   finite  difference  operators  represent- 

2  2 

9  9 

ing   -  — 2  ^^<^   ~  — 2   ^^  have 
3x         9y 

2 
(P+Q  V  ^,  =  M  Pv 
^   n+1        n 


or 


where 


P^/2,     .  ,,2  ^  pl/2  ^ 
n+1  n 


K  =  pl/2(P+Q)-^  pl/2 
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Thus 


II  P-"-/^    V       Jl     <    M^llKll     II  P-'-/^    V    II 
n+1      —  n 


and  since  K  is  Hermitian 


Kll  =  X    (K) 
max 


x,Kx 

max  -, — r 

(x,  x) 


(v,Py) 
max  ^^-^ — ^ — 


(y,Py)  +  (y,Qy) 


where 


Thus 


pl/2  y  =  K^/2  ^ 


P-"-/^  V   Jl  <  M^llp-"-/^  V  II 
n+1   —  n 


ow 


This  estimate  serves  to  indicate  that  for  subsonic  flows 
the  scheme  should  converge  at  a  rate  independent  of  the 
mesh  size. 

The  above  analysis  also  suggests  that  it  is  doubtful 
whether  such  a  scheme  would  converge  for  a  supersonic  fl 
with   M  >  1.   The  argument  presented  in  the  Appendix 
in  fact  indicates  that  the  scheme  would  definitely  not 
converge  for  a  linearized  supersonic  flow.   It  thus 
appears  that  the  fast  elliptic  solver  used  on  its  own  is 
not  likely  to  lead  to  good  conver  ence  when  the  supersonic 
zone  is  large.   If,  however,  it  could  be  supplemented  with 
another  scheme  which  give  fast  convergence   in  the  super- 
sonic zone,  the  two  in   combination  might  produce  an 
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effective  iterative  scheme.   In  fact  the  standard  line 
relaxation  method  for  transonic  flow  calculations  is  such 
a  scheme.   For  the  small  disturbance  equation,  or  in  the 
case  of  the  full  potential  equation  with  the  flow  aligned  with 
one  coordinate  direction,      the  method  consists  of  freez- 
ing the  nonlinear  coefficients  at  values  determined  from 
the  previous  iteration  and  solving  the  resulting  wave 
equation    in  the  supersonic  zone  by  a  marching  procedure. 
Thus  an  exact  solution  of  the  supersonic  zone  could  be 
obtained  in  one  step  if  the  correct  coefficients  and  data 
at  the  sonic  line  could  be  inserted. 

Thus  the  following  scheme  is  proposed:   use  a  two 
stage  iteration,  in  which  the  first  stage  is  a  step  using 
the  Laplacian  on  the  left-hand  side,  and  the  second  stage 
consists  of  a  fixed  number  of  relaxation  steps  to  stabilize 
the  supersonic  zone.   For  the  case  of  flow  aligned  with 
the  coordinate  system  a  single  relaxation  step  should  be 
sufficient.   The  full  potential  equation  in  a  curvilinear 
coordinate   system  with  an  arbitrary  flow  direction  requires 
the  use  of  a  rotated  upwind  difference  scheme  [3].   Then  a 
simple  marching  scheme  can  no  longer  be  used  in  the  super- 
sonic zone,  and  several  relaxation  steps  may  be  required 
in  the  second  stage. 
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2 .    Application  to  the  Transonic  Potential  Flow  Equation  in 

a  Mapped  Domain 

The  use  of  a  fast  Poisson  solver  requires  a  simple  domain 
such  as  a  rectangle.   This  leads  to  a  difficulty  in  applying 
the  proposed  method  to  an  exterior  flow  problem  with  an  infinite 
domain.   This  can  be  circumvented  by  using  the  full  potential 
equation  and  mapping  the  exterior  of  the  profile  onto  the 
interior  of  a  circle.   If  2ttE  is  the  circulation  it  is  conveni- 
ent to  use  a  reduced  potential  G  defined  by 

4,  =  G  +  ^°^-i-  E  6. 
^         r 

Then  G  is  finite  and  single  valued.   Now  a  fast  solver  for 
Poisson 's  equation  in  polar  coordinates  can  be  used  in  the 
first  stage  of  the  iteration.   For  this  purpose  a  scheme  using 
the  Buneman  algorithm  in  the  9  direction  has  been  programmed. 

Two  variants  of  this  approach  have  been  tried.   The  first 
treats  the  potential  equation  in  quasilinear  form.  The  residual 
at  each  point  is  evaluated  as 

R  =  (a^-v^)  r  1^  (r  G^)  +  (a^-u^)  Gq  g  -  2uv  r(Gj.Q+Gg-  E) 

+  (u^-v^)r  G^  +  (u^+v^)  (^  H,  +  V  H^)  , 
r  r   fc)      r 

where  H  is  the  modulus  of  the  transformation  to  the  exterior 

of  the  circle,  u  and  v  are  the  velocity  components 

2 

r(GQ-E)     -    sin    0  r   G      -    cos    9 

u  = ,         V  = 


H  >  -  jj 

and      a      is    the    speed   of   sound.      If      y      is    the    ratio   of 
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specific  heats,   a   is  determined  from  the  stagnation  speed 
of  sound   aQ   by  the  relation 

2     2    j^    ,    2^  2. 
a  =  a-j  -  -^-2—  (u  +v  )  . 

In  evaluating   R   upwind  differencing  is  used  in  the  usual 
manner  at  supersonic  points.   In  the  first  stage  of  the 
iteration  the  correction   C   is  determined   by  solving 

-  h   (-  ^r^  ^  Se  =  \ 

3. 

and   then   G   is    updated  by    the    rule 

G      =    G   +    CO    C 

where  the  superscript   +   denotes  the  new  value,  and 

0)   is  an  overrelaxation   factor.   In  the  second  stage  of 

the  iteration  an  ordinary  relaxation  step  is  used. 

Results  with  this  approach  have  been  quite  promising. 
Numerical  tests  have  confirmed  that  the  scheme  sometimes 
diverges  when  the  relaxation  step  is  not  included.  When  it  is 
included   fast  convergence  has  been  obtained.   Figures  1 
and  2   show  typical  results.   In  each  calculation  the 
calculation  was  performed  first  on  a  mesh  with  64  cells  in 
the  9  direction  and   16  cells  in  the  r  direction,  and  then 
on  a  mesh  with  12  8  x  32  cells.   The  interpolated  coarse  mesh 
solution  was  used  as  the  starting  point  for  the  fine   mesh 
calculation.   The  largest  absolute  value  of  the  residual 
anywhere   in  the  field  was  used  as  a  measure  of  convergence. 
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The  first  example  is  the  64A410  airfoil  at  Mach  .720. 

-1       -9 
In  this  case  the  residual  was  reduced  from  '^^10    to  10 

-3      -9 
in  26  cycles  on  the  coarse  mesh,  and  then  from  '^lO    to  10 

in  21  cycles  on  the  fine  mesh,  each  cycle  consisting  of 

one  Poisson  step  plus  one  relaxation  step.   The  Poisson 

step  takes  about  the  same  time  as  2  relaxation  steps,  so 

each  complete  cycle  requires  about  the  same  time  as  3 

relaxation  steps.   On  the  CDC  6600  at  the  ERDA  Computing 

Facility  at  New  York  University  one  complete  cycle  on  the 

fine  mesh  takes  about  1.5  seconds.   The  entire  calculation 

for  the  6  4A410  took  48  seconds.   The  second  example  shows 

a  shock-free  supercritical  airfoil  designed  by  Garabedian  [6] 

In  this  case  27  cycles  were  required  to  reduce  the  largest 

-9 
residual  to  10    on  the  coarse  mesh,  and  another  2  7  cycles 

-9 
to  reduce  it  to  10    on  the  fine  mesh.   In  corresponding 

calculations  using  relaxation  steps  without  the  Poisson 

steps  the  largest  residual  was  still  '^^10    after  1000  cycles 

on  the  fine  mesh. 

The  second  variant  treats  the  potential  equation  in 

conservation  form  using  a  rotated  difference  scheme  in  the 

supersonic  zone  [5].   In  this  case  the  residual  is  evaluated 

as 

R  =  r  1^  (rpV)  +  1^  (PU) 


where 


U  =  G,  -  E  -  ^^    ,    V  =  r  G^  -  ^^  , 


and   p   is  the  density.   If   M   is  the  free  stream  Mach 
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number   p   is  determined  from  the  speed  of  sound   a 
by  the  relation 

p  '    =  M   a   . 

^  oo 

Now  in  the  first  stage  of  the  iteration  the  correction  C 
is  determined  by  solving 

r  |-  (r  C  )  +  C„„  =  -  . 
9r     r      66    p 

The  second  stage  consists  of  k  relaxation  cycles.  Typically 
k  =  5. 

In  this  case  also  the  combined  iteration  has  proved 
to  give  faster  convergence  than  the  simple  relaxation 
method.       The  improvement  is  not  as  great  as  with  the 
first  variant,  however,  because  of  the  need  to  use  more 
relaxation  steps  than  Poisson  steps.   As  an  example  of 
the  application  of  the  method.  Figure  3  shows  the  pressure 
distribution  for  the  6  4A410  at  Mach  .720  recalculated  using 
conservation  form.    The  proper  theoretical  jump  condition 
is  now  satisfied,  as  can  be  seen.   In  this  calculation 

the  number  of  cycles  required  to  reduce  the  largest  residual 

-9 
to  10    was  37  on  the  coarse  mesh  and  40  on  the  fine  mesh. 

Each  cycle  consisted  of  1  Poisson  step  plus  5  relaxation 
steps,  and  took  about  the  same  amount  of  time  as  7  relaxa- 
tion steps ,  so  the  fine  mesh  calculation  is  equivalent  to  a 

little  under  300  relaxation  steps,  which  would  be  enough 

-5 
to  reduce  the  largest  residual  to  '^'lO    using  relaxation 

alone. 


3.    Conclusion 

It  is  concluded  that  the  use  of  a  fast  elliptic  solver 
in  combination  with  relaxation   is  an  effective  way  to 
accelerate  the  convergence  of  transonic  flow  calculations, 
particularly  when  a  marching  scheme  can  be  used  to  treat 
the  supersonic  zone  in  the  relaxation  process.   Other 
methods  of  preventing  divergence  in  the  supersonic  zone 
should  be  investigated.   Possibly  it  would  be  sufficient 
to  sweep  only  the  supersonic  zone   in  the  second  stage  of 
the  iteration.   This  would  lead  to  further  useful  savings 
of  computer  time. 
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Appendix.   Analysis  of  the  Poisson  Iteration  for  the 
Linearized  Equation  with   M  >  1. 

Let  the  equation 

(1-M^)  (J)    +   <i>        =  0  , 

2 

with  M   a  constant  >  1 ,  be  approximated  with  equal  mesh 

spacing  in  the  x  and  y  directions  by  the  Murman  difference 
scheme 


(1-M  )  ((}).  .-2<t>.     .     .  +  (J).  -  .)  +  4).  .^    -    2<p.  .    +    (p  .      ^    =  0 
ID    i-l/D   1-2, J      i,D  +  l      ID     I'D-l 


in  which  an  upwind  difference  formula  is  used  for  (j) 


XX 


Denoting    updated   values   by    the    superscript   +,    consider   the 
iteration 


4) 


+  +  +  + 

1+1, J  '^i]  1-1, D  i/D 


1+1  *1D  i/D-1 


=    <p  .  ^.      .    -    2^.  .    +    <^.     ,      .+     (M-l)((}i..-2(}).     ,      .  +  (ti.     „     .). 
^1+1,3  ^ID  i-l/D  ID  i-l/D       1-2, D 

Let   i    and    j    both    run    from   1    to  n    and   define    the   n    x   n    matrix 


T    = 


■2 

1 


1 
■2 

1 


1 
-2         1 
1-2         1 


Also  define  the  n  ^  n  matrix 
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R      = 


1 
■2 

1 


1 
-2 

1 


1 
-2 


and    let      $      be    the    matrix  with   entries      'J';-;- 
iteration   can   be   written    as 


where 


T$      +(J>T=(T+    aR)     $ 


a   =   M    -1    >    0 


Then    the 


Also   it   is   easily   verified   that 


where 


S      = 


0 

1 


R   =    T    S 


0 

1         0 

1         0 


-l/(n+l) 
-2/(n+l) 
-3/(n+l) 
-4/(n+l) 


Thus 


If 


T    $"*"    +    ^'^T    =   T  (I    +    aS)     <l> 


$      =    X    $ 


then      f    (with   its   elements    suitably   ordered   to    form   a   vector) 
is    an   eigenvector   of   the    iteration  matrix.      Consider   the 
form 

$   =    uv 
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where   v   is    an   eigenvector   of   T 

T   V  =    y   V 
Then      $      is    an   eigenvector   if 


T(I    +    aS)uv         =    A(Tuv      +    uv    T) 


=    A(T   +    yl)     uv'^ 


This  is  satisfied  if 


where 


P^u  =  X  u 


P   =  (T+pI)    T(I+aS) 


Thus  the  eigenvectors  of  the  iteration  matrix  can  be 

expressed  as 

T 
u  V 

where   v  is  an  eigenvector  of  T  with  eigenvalue   p   and 

u  is  an  eigenvector  of   P   ,  and  the  eigenvalues  of  the 

iteration  matrix  are  the  eigenvalues  of   P    for 

y  =  y,,y2/.-./y  •   For  large   n   the  smallest  eigenvalue 

y  of  T  is  of  order  —j-      But  as   y  ->■  0   the  eigenvalues 

n 

of  P   approach  those  of 

y 

P   =  I  +  aS 

and  correspondingly  some  of  the  eigenvalues  of  the  itera- 
tion matrix  approach 

1  +  a  A. 

1 


where  A.  are  the  eigenvalues  of   S, 


Now 
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J  1.    /  ^ -r  ,-.\    i^  .   n   ,n-l      ,   1 
det  (XI-S)  =    X      +   —-rr   X  ...    + 


n+1 n+1 

Thus  if  the  polynomial 

J   =  (n+1)  a"  +  n  X""-"-   .  .  .  +  1 
n 

has  a  root  in  the  right  half  plane   the  corresponding 
eigenvalues  of   Pq   will  lie  outside  the  unit  circle. 
Applying  the  Routh  Hurwitz  test,   it  can  be  verified 
that   J    has  at  least  one  root  in  the  right  half  plane 
when  n  >  4 . 


-17- 


This  report  was  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Administration, 
nor  any  person  acting  on  behalf  of  the 
Administration: 

A.  Makes  any  warranty  or  representation, 
express  or  implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  Information  contained  in  this  report, 
or  that  the  use  of  any  information, 
apparatus,  method,  or  process  disclosed 
in  this  report  may  not  infringe  privately 
owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  information,  apparatus, 
method,  or  process  disclosed  in  this 
report . 

As  used  in  the  above,  "person  acting  on  behalf 
of  the  Administration"  includes  any  employee 
or  contractor  of  the  Administration,  or 
employee  of  such  contractor,  to  the  extent 
that  such  employee  or  contractor  of  the 
Administration,  or  employee  of  such  contractor 
prepares,  disseminates,  or  provides  access  to, 
any  information  pursuant  to  his  employment  or 
contract  with  the  Administration,  or  his 
employment  with  such  contractor. 


-18- 


